\section{统计研究与随机模拟}

现代统计学期刊发表的论文中一半以上包括随机模拟结果 (也叫数值结果)。随机模拟可以辅助说明新模型或新方法的有效性, 尤其是在很难获得关于有效性的理论结果的情况下让我们也能说明自己的新方法的优点。

{\color{red}\textbf{[关键理解]:} 随机模拟在统计研究中扮演``虚拟实验室''的角色——当数学推导难以完成时,我们可以通过大量重复试验来探索方法的性能,就像物理学家用实验验证理论一样。}

比如, 为了掌握统计量的抽样分布经常需要依靠随机模拟。在独立正态样本情况下, 我们已经有样本均值和样本方差的分布,但是对其它分布 \(F\left( x\right)\) 的样本,我们一般只能得到统计量 \(\widehat{\theta }\) 的大样本性质,在中小样本情况下很难得到 \(\widehat{\theta }\) 的抽样分布理论结果,随机模拟可以解决这样的问题。在抽取 \(F\left( x\right)\) 的很多批样本后,得到很多个统计量 \(\widehat{\theta }\) 样本值,从这些样本值可以估计 \(\widehat{\theta }\) 的抽样分布,并研究样本量多大时 \(\widehat{\theta }\) 的抽样分布与大样本极限分布相符,计算估计的均方误差, 等等。

{\color{red}\textbf{[实际应用]:} 这就像重复抛硬币——单次结果不可预测,但重复百万次就能精确估计概率分布。对于复杂的统计量,理论推导可能极其困难,但重复抽样\(M\)次就能用经验分布逼近真实分布。}

下面举一个置信区间比较的例子说明统计研究中随机模拟的做法。实际应用中经常需要计算某个百分比的置信区间, 比如北京市成年人口中希望房价降低的人的百分比的置信区间。

{\color{red}\textbf{[研究问题]:} 这个例子展示了随机模拟的典型用途——比较两种竞争方法的优劣。当理论分析难以判断哪种置信区间更好时,我们设计对照实验,让数据说话。}

设随机抽取了样本量为 \(n\) 样本,其中成功个数为 \(X\) 个,则 \(X \sim  \mathrm{B}\left( {n,p}\right) ,p\) 为成功概率,要计算 \(p\) 的置信度为 \(1 - \alpha\) 的置信区间。记样本中成功百分比为 \(\widehat{p} = X/n\) 。假设样本量 \(n \geq  {30}\) 。

当 \(n\) 较大时由中心极限定理可知

\[
W = \frac{\widehat{p} - p}{\sqrt{p\left( {1 - p}\right) }/\sqrt{n}} \tag{3.70}
\]

近似服从标准正态分布。因为 \(\widehat{p} \rightarrow  p\) a.s. \(\left( {n \rightarrow  \infty }\right)\) ,所以把(3.70)式分母中的 \(p\) 换成 \(\widehat{p}\) ,仍有

\[
Z = \frac{\widehat{p} - p}{\sqrt{\widehat{p}\left( {1 - \widehat{p}}\right) }/\sqrt{n}}
\]

近似服从标准正态分布。于是可以构造 \(p\) 的置信度 \(1 - \alpha\) 的置信区间为

\[
\widehat{p} \pm  \frac{\lambda }{\sqrt{n}}\sqrt{\widehat{p}\left( {1 - \widehat{p}}\right) } \tag{3.71}
\]

其中 \(\lambda  = {z}_{1 - \frac{\alpha }{2}}\) 是标准正态分布的 \(1 - \frac{\alpha }{2}\) 分位数。当 \(n\) 很大时这个置信区间应该是比较精确的,即 \(p\) 落入置信区间的实际概率应该很接近于标称的置信度 \(1 - \alpha\) ,但是在 \(n\) 不是太大的时候, \(p\) 落入置信区间的实际概率就可能与标称的置信度 \(1 - \alpha\) 有一定的差距,我们要用随机模拟考察 \(1 - \alpha ,n,p\) 的不同取值下置信区间(3.71)的覆盖率。称置信区间实际包含真实参数的概率为覆盖率。

参数 \(p\) 的另一个置信区间,称为 Wilson 置信区间,是直接求解关于 \(p\) 求解不等式

\[
\left| \frac{\widehat{p} - p}{\sqrt{p\left( {1 - p}\right) }/\sqrt{n}}\right|  \leq  \lambda  \tag{3.72}
\]

记

\[
\widetilde{p} = \frac{\widehat{p} + \frac{{\lambda }^{2}}{2n}}{1 + \frac{{\lambda }^{2}}{n}},\;\delta  = \frac{\lambda }{\sqrt{n}}\frac{\sqrt{\widehat{p}\left( {1 - \widehat{p}}\right)  + \frac{{\lambda }^{2}}{4n}}}{1 + \frac{{\lambda }^{2}}{n}},
\]

则求解(3.72)得到的 \(p\) 的置信区间为

\[
\widetilde{p} \pm  \delta  \tag{3.73}
\]

我们用随机模拟来考察(3.71)和(3.73)在不同 \(n,p\) 下的覆盖率并比较这两种置信区间。

好的置信区间应该满足如下两条要求:

(1)覆盖率大于等于置信度,两者越接近越好；

(2)置信区间越短越好。

下面设计这个模拟试验。设总共重复试验 \(M\) 次。 \(M\) 越大,对覆盖率和置信区间长度期望估计的随机误差越小。比如我们希望模拟计算的覆盖率误差控制在 0.1\% 以下 (在 95\% 置信度下)。对每组样本, 真值是否落入计算得到的置信区间是一个成败型试验, 设真实覆盖率为 \(r,M\) 次重复试验中置信区间包含真实参数的次数为 \(V\) ,则 \(V \sim  \mathrm{B}\left( {M,r}\right)\) ,覆盖率的估计值为 \(V/M\) ,此估计的标准误差为 \(\sqrt{r\left( {1 - r}\right) }/\sqrt{M}\) 。如果覆盖率 \(r\) 近似等于置信区间的置信度 \(1 - \alpha\) ,则若 \(1 - \alpha  \geq  {0.8}\) ,近似地有 \(r\left( {1 - r}\right)  \leq  {0.8} \times  {0.2} = {0.16},M\) 次试验的覆盖率估计的标准误差为 \(\sqrt{0.16}/\sqrt{M} = {0.4}/\sqrt{M}\) ,以 2 倍标准误差作为覆盖率估计的误差大小界限(根据中心极限定理可知有 \({95}\%\) 的概率使得 \(\left| {\frac{V}{M} - r}\right|\) 小于 2 倍标准误差),令 \(2 \times  {0.4}/\sqrt{M} \leq  {0.001}\) 则有 \(M = {640000}\) ,这个试验重复量很大,如果程序耗时过长我们只好降低对随机误差幅度的要求。对于更复杂的问题,如果难以得到所需重复次数 \(M\) 的理论值,可以逐步增加 \(M\) ,直到结果在需要的精度上不再变化为止。

{\color{red}\textbf{[量级估计]:} 为达到0.1\%的精度需要\(M=640000\)次重复——这看似巨大,但正是``大数定律''的力量:标准误差以\(\sqrt{M}\)下降,要将误差减半需4倍样本量。这就是为什么高精度模拟计算量惊人。}

比较 (3.71) 和 (3.73)要考虑多种 \(\left( {n,p,1 - \alpha }\right)\) 组合的影响。取 \(1 - \alpha  = {0.99},{0.95},{0.90},{0.80}\) 几种。(n, p)的值要考虑到各种不同情况,初步选取如下的一些(n, p)组合:

\[
n = {30},p = {0.1},{0.3},{0.5},{0.7},{0.9};
\]

\[
n = {120},p = {0.05},{0.1},{0.3},{0.5},{0.7},{0.9},{0.95}
\]

\[
n = {480},p = {0.01},{0.05},{0.1},{0.3},{0.5},{0.7},{0.9},{0.95},{0.99}
\]

这样,我们设计了 \(\left( {n,p,1 - \alpha }\right)\) 的 \(4 \times  \left( {5 + 7 + 9}\right)  = {84}\) 种不同情况,每种可以得到四个数值:置信区间(3.71)是否包含参数真值 (1 为包含, 0 为不包含), 置信区间(3.73)是否包含参数真值,置信区间(3.71)的长度,置信区间(3.73)的长度。对每种给定 \(\left( {n,p,1 - \alpha }\right)\) 的情况,我们重复试验 \(M\) 次,分别得到置信区间 (3.71) 的覆盖率 \({\widehat{r}}_{1}\) ,置信区间 (3.73) 的覆盖率 \({\widehat{r}}_{2}\) ,置信区间(3.71)的平均长度 \({\bar{l}}_{1}\) 和长度标准差 \({s}_{1}\) ,置信区间(3.73)的平均长度 \({\bar{l}}_{2}\) 和长度标准差 \({s}_{2}\) 。

{\color{red}\textbf{[实验设计]:} 84种情况的全面测试体现了统计研究的系统性——不能只测试``理想情况'',必须覆盖小样本\((n=30)\)、极端概率\((p=0.01, 0.99)\)等``边界条件'',因为实际应用中这些情况常常出现。}

实际编程时, 我们先编写 84 种不同情况中的一种情况去试验, 首先试验较少的重复数, 比如 \(M = {100}\) 。当程序检查无误后,再逐步增大重复次数看计算时间和存储能力是否可以接受。对此例发现 \(M = {640000}\) 可行。最后,用循环处理所有 84 种情况,每种情况重复试验 \(M = {640000}\) 次。最后的结果汇总成表3.1(限于篇幅,结果有删减)。表中每一行是一种 \(\left( {1 - \alpha ,n,p}\right)\) 的组合。注意, \({\bar{l}}_{1}\) 和 \({\bar{l}}_{2}\) 的精度可以用其标准误差 \({\mathrm{{SE}}}_{1} = {s}_{1}/\sqrt{M}\) 和 \({\mathrm{{SE}}}_{2} = {s}_{2}/\sqrt{M}\) 来估计。

从表3.1的结果看出,即使样本量已经很大 \(\left( {n = {480}}\right)\) ,公式 (3.71) 的覆盖率仍然会低于置信度 \(1 - \alpha\) ,特别是当 \(p\) 靠近 0 或 1 的情况下公式(3.71)给出的置信区间更差。另一方面,公式(3.73)的覆盖率除少数例外总是略高于标称的置信度而且与标称置信度差距一般不大,两个公式得到的置信区间的平均长度相近。但是, \(1 - \alpha  = {0.8},n = {480},p = {0.01}\) 的 Wilson 置信区间覆盖率只有 74.7\%,说明 Wilson 置信区间仍有改进必要。

{\color{red}\textbf{[模拟揭示的真相]:} 这个结果令人警醒——即使\(n=480\)的``大样本'',传统方法(3.71)在极端情况下仍表现不佳,覆盖率系统性偏低。若仅靠数学推导,这些问题很难被发现;正是大规模模拟暴露了理论与实际的差距。}
